/*****************************************
separse.c

	spin echo parse sequence
	adapted from palm.c from Rose.c from episs.c
	waveforms are pre-computed externally in matlab or your favorite enviroment

	written 24 Feb 2009 -- msb and nzl

	*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
	sequence timeline.

	wait for TR - seqtime
	optional fat sat block - not added yet...
	excitation block
		turn on slice select gradient, gss
		wait for rise gss rise time
		wait rof1, pi/2 pulse transmitter for p1 ms, wait rof2
		turn off all gradients
		wait for gss fall time
		turn on refocus pulse, gssr
		wait for rgss rise
		wait for tssr
		turn off refocus gradient
		wait for fall time, rgssr
	acquire fid data
		while playing ro and pe gradients for te_delay time
			collect fid data
	spin echo block
		turn on slice select, and ro and pe crushers
		wait for slice select and crusher rise time
		wait for tcrush 
		turn off crusher
		wait for crusher fall time
		wait rof1, pi pulse for p2 ms, wait rof2
		turn on crusher
		wait for crusher rise
		wait for tcrush
		turn off all gradients
		wait for gradient fall
	acquire echo
		collect data for (at-te_delay)

	*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*/

#include "standard.h"
#include "group.h"
#define MAX_DAC_VAL 32767.0 /* system max DAC value */

	pulsesequence()
{
	/* our variables for this sequence */
    double  sefoo;                                      /* np fudge factor WTF? */
	double  seqtime;									/* sequence length */
	double  rgss;										/* ramp time of slice select gradient */
	double  tssr;										/* time of slice select refocus */
    double  tcrush1,tcrush2;                            /* crusher length */
	double	tcrushf;									/* crusher fraction, tweaker */
    double  rcrush;                                     /* crusher rise time */
	double  fid_dur,s_e_dur; 							/* length of the shaped gradients in seconds*/
	double	te_dur1,te_dur2,te_dur3;					/* the three parts of the te delay. */
    double  grate;                                      /* gradient slew rate */
    double  predelay;                                   /* delay before execution for TR */
	double		fid_np,s_e_np;								/* number of acq pts before and after 180 */
	char    ro_fid_file[MAXSTR], pe_fid_file[MAXSTR];	/* fid gradient table file path */
	char    ro_s_e_file[MAXSTR], pe_s_e_file[MAXSTR];   /* spinecho gradient table file path */
    char    tempstr[MAXSTR];
	char    *home_dir; 
	char	calib[MAXSTR];								/* calibration flag */
    

	/* bring in varian defined vars */
	initparms_sis();
	/* select varian variables ***********************
	tr		repetition time of sequence
	te		echo time in seconds
	gmax	maximum gradient strength
	trise	rise time to gmax
	gss		slice select gradient strength
	np		number of points of data to acquire
	gpe		phase encode gradient str
	gro		read out grad str, should equal gpe
    tcrush  crusher length
	**************************************************/

	/* load the gradient tables */
	if ((home_dir = (char *) getenv("HOME")) == NULL) 
	{
		abort_message("Error in determining HOME environment variable\n");
	}
	getstr("ro_fid_file",tempstr);
    sprintf(ro_fid_file, "%s.GRD",tempstr);
	getstr("pe_fid_file",tempstr);
    sprintf(pe_fid_file, "%s.GRD",tempstr);
	getstr("ro_s_e_file",tempstr); 
    sprintf(ro_s_e_file, "%s.GRD",tempstr);
	getstr("pe_s_e_file",tempstr);
    sprintf(pe_s_e_file, "%s.GRD",tempstr);



	/* get any other variables from the environment */
	fid_dur = getval("fid_dur"); /* get the shaped gradient durations */
	s_e_dur = getval("s_e_dur");
	tcrushf = getval("tcrushf"); /* get crusher tweaker */
    sefoo   = getval("sefoo");
    getstr("calib",   calib);    /* for slice selctive cal with big phantoms */ 
	
	/* calculate constants */
	grate   = trise/gmax;      /* gradient rate */
	rgss    = grate*fabs(gss); /* gss rise time */
	rcrush  = grate*fabs(gss); /* crush rise time */
	gcrush  = gss;             /* use gss as crush strength for now */
	tcrush1 = tcrush;          /* pre 180 crusher duration */
	tcrush2 = tcrush*tcrushf;  /* post 180 crusher duration */
	
	/* time and np bookkeeping ***********************************
	seqtime		total time needed to execute sequence
	tssr		time of slice refocus pulse plateau, set to make refocus integral = 1/2 slice select integral
	te_dur1		time from middle of excite to start of fid acq
	te_dur2		time of fid acq
	te_dur3		time from end of fid acq until middle of 180 pulse
	predelay	time to wait until execution, padding for TR
	fid_np		np before 180
	s_e_np		np after 180
	*/
    tssr     =  p1/2.0 + rof2 - rgss/2.0;
	seqtime  = 	6.0*rgss + 2.0*rof1 + p1 + 2.0*rof2 + tssr + at + p2 + tcrush1 + tcrush2 + 2.0*rcrush;
	te_dur1  =  p1/2.0 + rof2 + 3.0*rgss + tssr;
	te_dur3  =  rgss + rof1 + p2/2 + tcrush + rcrush;
	te_dur2  =  1000.0*te/2.0 - ( te_dur1 + te_dur3 );
	if (te_dur2 < 0.0)
	{
		abort_message("%s: Requested te too short. Min te = %f ms",seqfil,te_dur1+te_dur3);
	}
	
	predelay =  tr - seqtime; 
	if (predelay < 0.0) /* is TR long enough? */
	{
		abort_message("%s: Requested tr too short. Min tr = %f ms",seqfil,seqtime*1000.);
	}
	fid_np   = floor(te_dur2*sw*2.0/1000.0);
	s_e_np	 = np - fid_np + sefoo;
    
    /* diagnostics ***************************************/
    printf("****************\n");
    printf("sw:     %.3e\n", sw);
    printf("np:     %f\n",   np);
    printf("calib:  %s\n",   calib);
    printf("fid_np: %f\n",   fid_np);
    printf("s_e_np: %f\n",   s_e_np);
    printf("fid_dur:%f\n",   fid_dur);
    printf("s_e_dur:%f\n",   s_e_dur);
    printf("grads:  %s\n",   ro_fid_file);
    printf("grads:  %s\n",   pe_fid_file);
    printf("grads:  %s\n",   ro_s_e_file);
    printf("grads:  %s\n",   pe_s_e_file);
    
	/* PULSE SEQUENCE ************************************
	******************************************************/ 
	status(A);
	/* Phase cycle: Part I *******************************/
	mod2(ct,v3);                    /*0101 0101*/
	dbl(v3,v3);                     /*0202 0202*/
	hlv(ct,v4);                     /*0011 2233*/
	mod2(v4,v4);                    /*0011 0011*/
	add(v3,v4,v1);                  /*0213 0213*/
	assign(v1,oph);
	add(one,v1,v2);
	assign(zero,v10);
	assign(zero,v12);
	/* tr delay **************************************************/  
	delay(predelay);       
	msloop(seqcon[1],ns,v13,v14);  /* Slice loop **********/       
	rotate();                      /* hardware based coordinate frame rotation */
	poffset_list(pss,gss,ns,v14);  /* set slice offset frequency */	
	/* 90 degree RF excitation section ***************************/ 
	observepower(tpwr1);
	if (calib[0] == 'r') obl_gradient(gss,0.0,0.0);
	if (calib[0] == 'p') obl_gradient(0.0,gss,0.0);
	if (calib[0] == 'n') obl_gradient(0.0,0.0,gss);
	delay(rgss);
	shapedpulse(p1pat,p1,v1,rof1,rof2); /* actual excitation */
	zero_all_gradients();
	delay(rgss);	
	/* slice refocus section *************************************/
	if (calib[0] == 'r') obl_gradient(-gss,0.0,0.0);
	if (calib[0] == 'p') obl_gradient(0.0,-gss,0.0);
	if (calib[0] == 'n') obl_gradient(0.0,0.0,-gss);
	delay(rgss);
	delay(tssr);                                 /* refocus time */
	zero_all_gradients();
	delay(rgss);
	obsoffset(resto);                            /* put the reciever back on freq */
	/* fid acquisition section *************************************/
	if (calib[0] == 'r') obl_shapedgradient(ro_fid_file,"","",fid_dur,gro,0.0,0.0,1,NOWAIT);
	if (calib[0] == 'p') obl_shapedgradient("",pe_fid_file,"",fid_dur,0.0,gpe,0.0,1,NOWAIT);
	if (calib[0] == 'n') obl_shapedgradient(ro_fid_file,pe_fid_file,"",fid_dur,gro,gpe,0.0,1,NOWAIT);
	acquire(fid_np,1./sw); /* acquire data from fid before spin echo */
	/* spin echo 180 degree RF section *****************************/
	obspower(tpwr2);
	poffset_list(pss,gss,ns,v14);                /* move transmitter to slice offset frequency */
	if (calib[0] == 'r') obl_gradient(0.0,gcrush,gcrush);
	if (calib[0] == 'p') obl_gradient(gcrush,0.0,gcrush);
	if (calib[0] == 'n') obl_gradient(gcrush,gcrush,0.0);
	delay(rcrush);
	delay(tcrush1);
	if (calib[0] == 'r') obl_gradient(gss,0.0,0.0);
	if (calib[0] == 'p') obl_gradient(0.0,gss,0.0);
	if (calib[0] == 'n') obl_gradient(0.0,0.0,gss);
	delay(rgss);
	shapedpulse(p2pat,p2,v2,rof1,rof2);          /* 180 excitation */
	if (calib[0] == 'r') obl_gradient(0.0,gcrush,gcrush);
	if (calib[0] == 'p') obl_gradient(gcrush,0.0,gcrush);
	if (calib[0] == 'n') obl_gradient(gcrush,gcrush,0.0);
	delay(rcrush);
	delay(tcrush2);
	zero_all_gradients();
	delay(rgss);
	/* spin echo acquisition section ********************************/
	obsoffset(resto);                           /* put the reciever back */
	if (calib[0] == 'r') obl_shapedgradient(ro_s_e_file,"","",s_e_dur,gro,0.0,0.0,1,NOWAIT);
	if (calib[0] == 'p') obl_shapedgradient("",pe_s_e_file,"",s_e_dur,0.0,gpe,0.0,1,NOWAIT);
	if (calib[0] == 'n') obl_shapedgradient(ro_s_e_file,pe_s_e_file,"",s_e_dur,gro,gpe,0.0,1,NOWAIT);
	acquire(s_e_np,1./sw);                     /* acquire data from spin echo */
	zero_all_gradients();
	endmsloop(seqcon[1],v14);
}

